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Abstract 

We present a lattice field theory of spins coupled to Dirac fermions, as a 
model for the doped copper oxide compounds. Both the fermionic and spin 
degrees of freedom are treated dynamically. The influence of the charge carri- 
ers on the magnetic ordering follows automatically. The magnetic phase dia- 
gram at zero temperature is studied numerically and with Mean-Field meth- 
ods. The Hybrid Monte Carlo algorithm is adapted to the 0(3) non-linear 
a model constraint. The charged excitations in the various phases are stud- 
ied at the Mean-Field level. Bound states of two charged fermions are found 
in a strongly coupled paramagnetic phase, without requiring a Fermi sea. 
We acquire a qualitative understanding of high-T c superconductivity through 
Bose-Einstein condensation of those dynamically bound pairs. The model 
also implies insulating behaviour at low doping, and Fermi liquid behaviour 
at large doping. We predict the possibility of reentrant superconductivity and 
the absence of superconductivity in spin-1 and spin-3/2 materials. 
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1 Introduction 



Since the discovery of the Cu-based high-T c superconductors Q there has been a lot of 
interest in the physics of antiferromagnets doped with holes. The understanding of such 
a system would allow a unifying picture of the magnetic and superconducting properties 
of those materials. Some of the properties to be considered are: 

1. Superconductivity appears in doped, ceramic materials, like La2Cu04 doped to 
La2- x -Sr x Cu04, in a narrow range of doping not too far from the antiferromagnetic 
insulating state, at low temperatures. 

2. For zero doping fraction, x = 0, they are insulating antiferromagnets. 

3. At large x they appear to have normal metallic behaviour. 

4. They are layered compounds, made up of Cu02 planes separated by some "charge 
reservoir" . The doping suppresses the antiferromagnetic (AFM) correlation between 
neighboring CuC>2 planes. However, large (£ ~ 10a) in-plane correlations remain 
just above T c . Also, transport phenomena occur mainly in the CuC>2 planes, so 
everything looks like a d=2 problem, with localized spins on the Cu and mobile holes 
on the O ions. 

5. As in the conventional superconductors, the superconductivity is charge-2 in nature, 
suggesting a pairing state. 

6. The normal {i.e. non superconducting) state is anomalous: there is experimental 
evidence for anomalous behaviour of several response functions (magnetic suscepti- 
bility, specific-heat, DC conductivity) in the normal phase. This is usually refered 
to as a pseudo-gap phase, for which a picture of a singlet pairing already above T c 
is emerging [||. That is, the pair formation seems to occur independently of the 
quantum liquid condensation (i.e. the macroscopic phase coherence), at least for 
underdoped (i.e. superconducting samples less than optimally doped) materials Q. 

7. There is a negligibly small isotope effect and there is evidence § for a d x i_ y i wave 
pairing state. 

The essential questions posed by the phenomenology of the new perovskite supercon- 
ductors have been phrased as follows in Ref. ||: 

• What is the physical origin of the anomalous normal state? 

• How can it be characterized? 

• What is the mechanism for high-T c superconductivity? 

• What is the pairing state? 

We will try to formulate an answer to these questions in the discussion section of this 
paper. 

The behaviour of the undoped compounds (point 2. above) is well understood in terms 
of the two-dimensional nearest-neighbour quantum antiferromagnetic (AFM) Heisenberg 
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model. This model can be mapped to the 0(3) non-linear cr-model ||, and the predictions, 
obtained through renormalization group analysis jjj and chiral perturbation theory |8), 
agree very well both with experimental |]] data for the cuprates and with numerical [jl(]] 
data from direct simulation of the AFM Heisenberg model. Therefore, it seems natural 
to take either the AFM Heisenberg model or the 0(3) sigma model as a starting point 
to study the doped materials. In contrast with the successful description of the undoped 
compounds, however, to our knowledge no model has yet been capable of describing the 
full doping range. 

We have recently proposed a simple model [11, 12 1 which is a natural extension of 
the 0(3) model to include doped charge carriers. The model is able to explain, at least 
qualitatively, many of the observed properties of the perovskites from undoped to (highly) 
overdoped compounds. It is a lattice-regularized, field-theoretical model of interacting 
spins and Dirac fermions in 2+1 dimensions, with only two free parameters in addition to 
the temperature: a nearest-neighbour spin coupling and a spin-fermion coupling. 

Conceptually, the model has much in common with microscopic spin-fermion Hamil- 
tonians which are somewhat less restrictive than the t-J model [14|. The formulation 
of our model as a local lattice field theory has several advantages over those and other 
models, though. Mean-field (MF) calculations are feasible and, in the absence of a sign 
problem, numerical Monte Carlo (MC) simulations can be carried out. Our MC results 
furthermore demonstrate that the MF approximation is quite reliable here. 

In the present article we want to present a careful, detailed discussion of the model, its 
symmetries, and its properties, and give full technical details and results of the MF and 
MC calculations, some of which were reported in Ref. JIl]]. In the present paper, our mean- 
field and numerical studies will be limited to the zero-temperature case, corresponding to 
infinite Euclidean time direction. However, we will argue on general grounds what happens 
when the temperature is increased. 

The remainder of this paper is laid out as follows. In Section 2 we present our model. 
In Section |2.1| we discuss the choice of lattice fermions, comment on the symmetries of the 
model, give its phase diagram and prove the reality of the fermion determinant, even in the 
presence of a chemical potential. In Section 2.2, we give our interpretation of the model 
as an effective field theory which tries to embody the essential features of holes strongly 
coupled to a dynamical AFM spin background, the only fundamental ingredients being 
symmetries. In Section 3 we examine the phase diagram of the model in the MF approxi- 
mation. In Section 4 we use MC simulations to complete the study of the phase diagram. 
For this purpose we have developed a new method that exactly solves the technical prob- 
lem related to the length- 1 constraint on the spin variables p5| . Section 5 is devoted to a 
study of the relevant excitations in the different phases of the system, at the MF level. A 
crucial result is the dynamical generation of spin singlet bosonic bound states of charged 
fermions in the so-called paramagnetic strong (PMS) phase. At the MF level we have 
not detected any light fermionic excitations at zero temperature in this PMS phase. An- 
other interesting result is the emergence of fermionic excitations around spatial momenta 
(±7r/2, ±vr/2) |l(| [l7]], in the low doping AFM phase preceding the superconducting (SC) 
phase. An important consequence of the paramagnetic nature of our SC phase is that 
fermionic excitations, if any, will not center around (±7r/2, ±7r/2), while such a behaviour 
is typical of an AFM phase (see ref. [18] for some recent experimental results). Next, in 
Section 6.1, we analyze how our model describes the various phenomena in the cuprates at 
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zero temperature and in Section 6.2 we conjecture the qualitative behaviour of the model 
at non zero temperature, using arguments based on a few very general thermodynamical 
properties. 

In the Discussion section (Section 7) we will try to answer the essential questions 
pointed out in Ref. H and quoted above. We also reconsider points 1 to 7 as well as the 
Bose-Einstein (BE) to BCS crossover from the perspective of our model, and we comment 
on our predictions. 

2 The model 

2.1 Motivation, Formulation, Symmetries, Phase Diagram 

In order to motivate our model, we start by collecting some of the most compelling pieces of 
experimental and theoretical evidence and turning them into guidelines for the formulation 
of the model. 

First of all, recall that the undoped parent compounds are described excellently by the 
2+1 dimensional 0(3) non- linear cr-model |7|, || [9|, (see also Ref. |l9|]). Therefore, it is 
very natural to take this as our first guideline: 

A. In the limit of zero doping, x = 0, our model should reduce to the 0(3) 
non-linear cr-model. 

Another important piece of information comes from point 7 of the list in Section 1. 
The virtual absence of an isotope effect and the evidence for a d x 2_ y 2 wave pairing state 
favor a mechanism based on spin rather than phonon-mediated pair formation |c], ||. 
This leads us to our second guideline: 

B: Our model should be a spin-fermion model describing fermions of varying 
mobility interacting with a dynamical background of localized spins. 

Next, as emphasized in |22|| , the general topology of the cuprates phase diagram was 
in fact first predicted in Ref. |23|] . Their point is to assume a strong coupling between the 
spin of the hole on the oxygen (see point 4) and the spin of the surrounding Cu 2+ ions 
due to the direct overlap of the wavefunctions. Therefore, from this: 

C: We expect the strong-coupling regime of our spin-fermion model to 
be the relevant regime to understand, eventually, the insulating and superconducting 
phases of cuprates. 

Notice that to in order to arrive at A and B we have essentially used only information 
from points 2, 4 and 7, while guideline C comes from the more general observation that 
we deal with a strong coupling phenomenon. 

Given the above considerations, we have proposed JiT| a model defined by the following 
(2 + l)-dimensional lattice euclidean (imaginary time) path integral, 

Z = J D<j) Dtp Dip exp(— 5) (1) 

with action, 

X,fl X,fl X 
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We use this expression as our starting point, but it should be noted that the model depends 
only on the ratio y = X/p, through a change in the normalization of the fermion field. In 
terms of the effective spin- fermion coupling y, we get: 

S = ~^2kcf) x - cf> x+fi + 2 V^7 M (?Ar+A - ipx-p) + 4>xcf> x ■ Tip x . (3) 

X,fl X,fl X 

Here x runs over a (2 + l)-dimensional cubic Euclidean space-time lattice. We keep in 
mind, however, that some coherence perpendicular to the CuC>2 planes is required to avoid 
problems with the Hohenberg-Mermin- Wagner theorem at finite temperature, and for a 
Bose-Einstein condensate to form. Some small supplementary couplings in the orthogonal 
spatial direction would play this role. 

The fields tp represent the doped electric charges. Each ip x is a fermionic four-spinor as 
a shorthand for two flavours of two-component Dirac spinors (we use flavour to mimic the 
spin, because there is no spin in two spatial dimensions). Each flavour accounts for two 
components of the four of a (3 + l)-dimensional Dirac electron or hole |24|] . Both flavours 
are taken in the same irreducible spinor representation, with 2x2 gamma matrices taken 
as the Pauli matrices a^. The 4x4 matrices 7^ in Eq. (^|) have the form 

^={ a o ^ = 1 > 2 ' 3 - & 

The kinetic term for the fermions is of the nearest-neighbour (hopping) form. Lattice 
fermions defined in this way undergo species doubling in the perturbative continuum 



limit [25]. But since all the fermions, the physical one as well as the doublers, decouple 
in the continuum limit performed in the region of strong spin- fermion coupling pf|, the 
doubling problem is in fact limited to the weakly coupled region. As our only relevant 
conclusion in this region, the occurrence of Fermi liquid behaviour because of the weak 
coupling, is not affected by the doublers, we canuse these "naive" fermions without any 
problem. The three-component fields 4> denote the spins located at the copper ions. They 
are real scalar bosonic variables, subject to the constraint (f) 2 = 1, as in the 0(3) non-linear 
cj-model. Their kinetic term, a nearest-neighbour hopping interaction, is the field-theoretic 
equivalent of a Heisenberg superexchange interaction. The last term in Eq. (|3|) describes 
the interaction between the spins and the Dirac fermions, which is diagonal in Dirac space. 
The Pauli matrices r a act in flavour space. 

Let us now consider the symmetries of (||). First of all, we have the usual cubic 
symmetry, the lattice remnant of (2+l)-dimensional Euclidean "Lorentz" symmetry. Next 
there is a discrete parity symmetry, which in (2+1) dimensions is defined as the reflection 
of one of the spatial axes, say the rc-axis. Under this parity symmetry, the fermions can 
be seen to transform as 

i>—KTiip, tp^-tpax, (5) 

so the 4> field is a pseudoscalar in this sense. In addition, the action (|||) is invariant under 
an SU(2) "flavour" symmetry in which i/j transforms as the fundamental representation 
and 4> transforms as the adjoint one. So, our model maintains after doping the well tested 
symmetry of the 0(3) non-linear cr-model describing the undoped material |J. Note that 
by requiring the two fermion flavours to have the same Lorentz structure (that is, by 
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choosing the 7's as in (Q)) no fermion mass term is allowed which preserves the above 
symmetries [27|. 



There are two more discrete symmetries of our model (||), which will be useful in the 
MF calculation of the phase-diagram. The first one is 

Z(k,y) = Z(k,-y), (6) 

which becomes clear if we make the change of variables 

ip x -> exp (i^ e xj i>%, $x exp (i^x) i>x, (7) 

where 

e n = (-l)^-M. (8) 

This implies that Z(k, y) is a function of y 2 only and we can restrict ourselves to y > 0. 
In addition, there is a symmetry 

Z(k,y) = Z(-k,-iy), (9) 

as can be seen by making the substitutions 

ip x -> exp (ij^j ipx, 4>x exp (ij^ i>x, 4> x ~* e x 4> x - ( 10 ) 

The latter symmetry implies that the lattice regularization of the non-linear cr-model, 
y=0 (or y=oo, see Sections 0, H), is equally valid in a ferromagnetic or an antiferromagnetic 



phase. This will be of fundamental importance in section |2.2 . 

In order to perform computations in field theoretical models of this type, one has to 
integrate out the fermion fields; the reason being that they are anticommuting Grassmann 
fields, while both in MF and in MC one needs to be able to work with a c-number or- 
dering. This integration leads to a ^-dependent determinant which is called the fermion 
determinant. It is important to know whether this determinant is a real number. To 
study this, let us write down the original fermion matrix (Latin letters x,y, . . . will refer 
to lattice points, will represent flavour indices, while Greek letters a, (3, ... are 

used for Dirac indices): 



M X ai;yf3j — ^-xai;y(3j ^xai;yPjj (H) 
o (5x+fi,y ~ Sx-n,y) Vg Sjj, (12) 



-^■xai;y/3j ~ ^ 

I" 

Y xa i;y/3j = V $xy {4> ' 1% ■ 6 a/3 . (13) 

Keeping in mind that for Pauli matrices 020^02 = — cr*, where * means complex conjuga- 
tion, and that [7, r] = 0, one easily proves that, for real y, 

a 2 r 2 (K + Y) a 2 T2 = - (K* + Y*) . (14) 

Therefore, 

det (K + Y) = det (-K* - Y*) = [det (K + Y)]* , (15) 
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i.e. the determinant is real. Thus, by doubling the number of fermion families, we obtain 
a positive determinant. Had we introduced a chemical potential, fj,, the only change would 
be the introduction of e ±/x on the temporal links of the kinetic matrix |2lJ . The essential 
requirement for Eq. (|l4|) to hold (that the only non-real numbers are in 7,r) is thus not 
endangered by the chemical potential and the determinant is still real. 




Figure 1: Phase diagram of the action (Q), for two fermion families. Dashed lines 
are from the MF calculation, solid lines from a MC calculation on an 8 3 lattice. The 
arrows indicate the doping- "trajectory". 

The phase diagram of the model at zero temperature is shown in fig. Notice that it 
is very similar to the phase diagram of (chiral) Yukawa models for the electroweak sector 



of the Standard Model of elementary particle interactions [29]. At y = oo and at y = we 
recover the non-linear a- model (see sections ^J4|) with its well known paramagnetic (PM), 
ferromagnetic (FM) and antiferromagnetic (AFM) phases. At finite y, we expect these 
phases to extend into the (k,y) plane. One of its most important features is that there are 
two mutually disconnected paramagnetic phases, one at weak coupling (called PMW) and 
one at strong coupling (PMS). One sees that the PMW-FM and the PMW- AFM transition 
lines meet in a point A, where this disordered phase ends. In the strong coupling sector of 
the phase diagram, a similar behaviour is found, with the two transition lines meeting at 
point B. This observation means that one may expect totally different behaviour in each 
of the two paramagnetic phases. This is indeed the case, as we shall see later. 
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As there is no evidence for a phase transition between the strong- and weak-coupling 
regions of the FM and AFM phases, we name them FM(W) and FM(S), AFM(W) and 
AFM(S) (note the parentheses). There may be crossovers between these regions, though. 

Between the points A and B, we find a phase where both the magnetization and the 
staggered magnetization are different from zero. We name this phase exotic magnetic 
(EM). An appealing possibility is that it corresponds to a helicoidal phase. This would 



merit a more detailed study (see section 3.1 for some comments on this point). We expect 



the EM phase to disappear for large enough —k, but we have not explored this numerically. 

2.2 Phenomenological use of the model 

Before embarking on a detailed study of the model, let us briefly discuss its relation with 
the phenomenology of the cuprates. In particular, it is important to explain how the 
doping fraction x is accounted for by the model. 

By increasing the doping fraction, one increases the overall mobility of the charge car- 
riers. In our model, this mobility is embodied in the fermionic hopping parameter p in 
(g). The exact mapping between x and p is not important for the qualitative description, 
but it is important that in the limit x, p — ► one recovers the 0(3) model of the undoped 



compounds (cf. guideline A. in Section 2.1). Thus, the undoped compound will be de- 



scribed by y = X/p = oo, and the large-?/ regime corresponds to the strongly-correlated 
small- x region, with immobile, localized carriers |30[| . 

Within the lattice 0(3) model, we still have two possible points to describe the undoped 
materials (see fig. |TJ) . One in the FM phase and one in the AFM phase, related by symmetry 
(0). Although one usually considers the FM phase, especially in continuum descriptions, 



it is perfectly legitimate (see Eqs. (|9|,10) and comments below it) to consider the AFM 
phase. In fact, as the doping is experimentally known to destroy the antiferromagnetism, 
this turns out to be the most natural option (see fig. |l|). (Let us remark though, that the 
reader can start in the ferromagnetic phase if he/she so wishes, by taking an imaginary 
coupling constant, as in Eq. @). 

The most important effect of changing the doping fraction, x, will be a change in the 
effective spin-fermion coupling, y (indicated by the arrow-line in figure |l|) , and therefore a 
change of the system ground state and excitations (fermionic and bosonic). Strictly speak- 
ing however, there will be a relation between x and a chemical potential, p, which in fact 
could be introduced in our formulation (recall that the fermion determinant would remain 
real). This relation would be in addition to the (phenomenological) relation between x 
and (y,k). We will now argue, though, that for the physics considered in this paper we 
can ignore the chemical potential. In other words, a possible chemical potential would not 
change the excitations of the ground state in any essential way. 

Let us start with a general point, that our model captures the insertion of charge 
carriers in an effective way, by increasing the overall mobility. One could object to this 
that a vanishing chemical potential entails the absence of a Fermi surface. In the strong- 
coupling region, however, there are two possible responses to this objection. First of all, 
for small to intermediate x, (i.e. before and just after superconductivity appears) the very 
concept of a Fermi surface is very problematic from the experimental point of view [|l^] . 
Second, it will turn out that the Fermi surface plays no role in our mechanism to form 
bound charge pairs. Therefore, ignoring a possible chemical potential in the small- x region, 
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whose main role would have been to fix the Fermi surface, is perfectly consistent, both 
with the experimental situation and within the context of the dynamical pair formation 
in our model. 

The overdoped, large- x region, on the other hand, is characterized by Fermi liquid 
behaviour, so the presence of a chemical potential in the small-y region of our model will 
be essential to determine the thermodynamics of the system. Since the only statement we 
want to make concerning this large- x region is the actual emergence of this Fermi liquid 
behaviour, we can forget about the chemical potential also in the large- x region, for our 
present purposes. 

In the superconducting regime, finally, a chemical potential would change to some 
extent the ground state and the excitations of the system. However, one expects that in 
a region of heavy fermions, a quenched region, this change will be small. So we can drop 
the chemical potential here too. 

After having discussed the phenomenological use and interpretation of the model de- 
fined by (||), we shall proceed to solve it. 



3 Mean Field Calculations of the Phase Diagram 

Our aim in this section is to determine the zero-temperature phase diagram of the model 
in the y-k plane (cf. Fig. ||), using Mean Field techniques. These calculations already 
provide a lot of insight, especially for the strong coupling region. They will be contrasted 
with numerical simulations for the phase diagram in Sect. ||, and they will be extended to 
a study of the relevant charged (quasi-particle) excitations in Sect. [|. 

Our Mean Field calculations are based on small- and large-y expansions combined with 
the saddle point methods described in Ref. ]3l|| . This approach guarantees a systematic 
expansion in 1/d, which is particularly important for operators which are zero to lowest- 
order. Our particular method furthermore allows us to handle (products of) fermionic 
variables occurring in the expansion of the fermion determinant in a well-defined way. 
These techniques have been developed and applied in the context of similar lattice models 



[[34 33] of the Electroweak sector of the Standard Model of elementary particle interactions. 

We shall first concentrate on the small-y region, and incorporate the fermion determi- 
nant up to 0(y 2 ). 

In order to apply the saddle-point method, the integration over the fields must be unre- 
stricted. We therefore need to replace the integration over the spin vectors cf>, constrained 
by the condition \cj)\ = 1, with an integration over unconstrained variables This is done 
by multiplying the functional integrand in Eq. (jl|) by 



n a=l J -°° 



nn 



roo j Aa 

/ ^exp[L42(#-g)] , 



x a=l ' 

and replacing a conveniently chosen subset of the <\> variables in the action S with £ fields. 
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We obtain 



D£DA 
3- exp 



x,/i a; 



(27T) 

DtpDip exp 



x,ii 



(16) 



Note that both the £ fields and the auxiliary fields A are unconstrained. 

Now we have to integrate out the constrained variables (f>® (as well as the fermions), 
before the mean fields can be introduced. Let us concentrate on a single 4> n integration 
in Eq. (|l6|), dropping the subscripts n for simplicity. First, we perform an expansion in 
powers of y. We can write 



— exp [iA ■ 4> - y 4> • rif] 



= exp [«(iA)] exp 
where we have defined 



-yQ° 



lA + 2y 2QaQbTab + °w 



(17) 



Q a = ^ T a ip , uUA) = In / ^ exp [iA ■ tf>] , T ab 

J Air 

and we have introduced the notation 



dcj> 



HAW hA> 



(0) iA = /^exp^.^ J%£ exp[iA.<t>}. 

In addition, we introduce a Hubbard-Stratonovich vector parameter A to deal with the 
quartic fermion term in Eq. (|17|), 



exp 



\y 2 £ Q a Q b T ab 



dX 



exp 



x exp 



(2tt) 3 / 2 



(18) 



ab 



(Note that the matrix T is self adjoint, and positive definite if A is imaginary, so the 
square root is well defined). Thus, up to this order in y 2 , the action is bilinear in the 
fermion fields. 

Carrying out the fermion integration in Eq. (|T6|) now gives det M, where 



ab 



M X)a ,r, yM = K Xta>i . yAj + yS xy 5 aP J2{(^)iA x - [V^J xb Kj ■ (19) 

a b 

The matrix K has been defined in Eq. (|l2[). 



10 



The mean fields are the field values at the saddle point of the free energy 



(20) 



x,ii 



A choice of the mean fields should be done at this point, as we cannot calculate logM for 
general {A x , X x }. An appropriate choice for the study of a PM-FM phase transition is 



in terms of which (N is the lattice volume) 



(0, 0, — ia) , 
(0,0, v), 
(0,0, A), 



(21) 



F/N 



-ma) 



kdv + av H — A 
2 



N 



Tr log M . 



with a, v and A satisfying the saddle point equations 



VF 



0. 



(22) 



(23) 



The fermion matrix, M(a,v, X), can be calculated in momentum space, where it is 
diagonal in its momentum indices. One easily finds 



det M = exp 



^ , El=i sm 2 p M + y 2 (u'(a) - Ay^R)' 



2 ^ 



E 



3 • 2 



(24) 



where we have divided out the determinant for free fermions. We need only the leading 
0(y 2 ) contribution to the exponent, hence the mean field free energy becomes, in the 
infinite volume limit: 



F/N 



where 



a 



(a) - kdv 2 + av + ^A 2 - 2y 2 (u'(a) - \^u"(a)^ C 
n d 3 p 



1 



1.0109240. 



(25) 



(26) 



Next, we shall discuss the actual solutions to Eqs. (p3[). From u{a) = In (sinh a/a), 
one easily finds that a = v = A = always fulfill them. For small k, y, it is a true 
minimum of the free energy. This characterizes a paramagnetic (PM) phase, since none 
of the fields develops an expectation value. 

For larger values of k and y, there is another, non-trivial solution, corresponding to a 
ferromagnetic (FM) phase. It emerges when a negative mode in F/N starts to develop, 
as a function of the mean fields, and the transition between the two regions is given by 
the condition (F" is the Hessian matrix) 



detF" 



( a =0,u=0,A=0) 



0. 



(27) 
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This condition is satisfied for F/N of Eq. (p5t) if 



This curve defines the phase transition line between the PM and FM phases in the 
small- y region. Using the symmetry (Q), we deduce that there is a similar transition 
separating the PM and AM phases, 

Let us finally remark that in the presence of Nf such fermion fields we would have Nf 
factors of det M, leading to a multiplication of Cq by Nf in Eqs. ( |28| ) and (p9|). 

The large-?/ region is easier to deal with. Here it is convenient to integrate out the 
fermions directly in Eq. ([16]), leading to (summation over repeated index is carried) 

det M^a^y^j = det ^K x>a ,i-,y } p,j + y5 a p8 xy ^ <fx T i^j (30) 
= det ^yS a ^S xz ^2^T^j 

x det [ 5 Z y8 lP 8 kj + $x T li K zakyfij ) ■ ( 31 ) 



y b 



Here we have used that (^2i a <fi a T a ) 2 = 1 (recall the </>'s are unit vectors). Now we can 
expand log(det M) in powers of 1/y. The 0{l/y) term vanishes by virtue of K xx = 0. To 
second order one obtains 

log det M = log y 4N + Tr £ ^*Wyi £ ^K tlj , y pA (32) 

= logj/^ + ^X^-^+A- (33) 

Here, log y 4N is an irrelevant constant that can be dropped. Notice also that this expression 
will acquire a prefactor Nf if there are Nf identical fermion flavours. One sees that, up 
to 0(l/y 2 ), the only effect of the fermion determinant is a renormalization of the scalar 
hopping parameter of the 0(3) model, 

k -» k + N f \ . (34) 

Note that we did not introduce any mean fields to derive this result. The usual MF 
treatment of the 0(3) model with this renormalized coupling now immediately gives us 
the required phase transition lines in the large-y region of our model: 
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It is interesting to compare the small- and large-y results, to leading order in 1/d. As 
is well known, the first order in this expansion is equivalent to any MF approximation, up 
to higher-order terms. For this purpose, we need the 1/d expansion of the constant Cq in 
Eq. (pq), which can be calculated as follows: 

CoW = r^ v , \ 2 = 2 r d s(e^I (s)r (36) 

-K i+ a +o (»))- 

where Iq(s) = J^ 7T (d6 /2tt) exp(scos#) is the modified Bessel function. In fact, the second 
equality in Eq. ([36]) was used to obtain the numerical result ( p6| ) for Co- 

Keeping only the leading-order term 2/d for Cq we find that the phase transition lines 
would meet at y 1 = 2/d. 

Now we are ready to map out the phase diagram of the model, as predicted by the MF 
method for the weak and strong coupling regions. This is done in Fig. |l[ The vertical axes 
at y = and y = oo correspond to the 0(3) model, with its disordered (PM) and ordered 
(FM and AFM) phases. These phases extend into the y-direction, both for y > and 
y < oo. Note that all the phase transition lines bend downward. This can be understood 
intuitively by assuming a MF value for the fermion condensate, which would act as an 
external field tending to align the spins 4> in parallel. 



4 Monte Carlo: Method and Results 



A well established method for dynamical fermion simulations is Hybrid Monte Carlo 
(HMC) |J4|. However, the implementation of this algorithm in a model with constrained 
variables is not straightforward. This has been satisfactorily achieved for models with vari- 
ables belonging to a Lie group [35], like SU(iV) gauge theories or like some spin-models, 
such as the 0(iV = 2,4) non-linear a-models. However, for other spin variables (not in a 
Lie group), as in the 0(3) non-linear cr-model, this had not been satisfactorily solved yet, 
although the problem arose already in the first simulations using the Langevin algorithm 



]. Our solution is a generalization of the strategy in [35]. 

We shall first discuss our solution in the quenched approximation, where comparison 

and then deal with the full theory in 



with other algorithms is possible (Section LI 



Section 4.2, Finally our Monte Carlo results for the phase diagram of the full theory will 



be presented in Section 4.3. 



4.1 The HMC method for the quenched approximation 

For the purpose of discussion it will prove convenient to briefly describe the HMC method 
for unconstrained bosonic variables 
presentation): 



(x) with action Sb(4>) (see ref. [36| for a pedagogical 
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1. Introduce uncorrelated gaussian variables n(x) of unit variance (the conjugate mo- 
menta for the fields 0) and define a Hamiltonian 

ff = £^r 2 (z)+S B (0). (38) 

X 

One can then use the hamiltonian equations of motion 

<fi(x,T) = 7t(x,t), (39) 
SS B 



7t(x, t) 



x,t) 



to perform a microcanonical Molecular Dynamics evolution in "Monte Carlo time" , 
t. After a certain period of MC time (called "trajectory"), new random momenta 
7r(x) are chosen ("refreshing" the momenta). The crucial properties of Eqs. ( f40|) are 
their time reversibility, and the invariance of the Liouville measure, Deft Dir, under 
time evolution. 

2. In practice, the molecular dynamics equations of motion for a trajectory are dis- 
cretized into N steps At. This is done using a leap-frog algorithm which is exactly 
time reversible, but does introduce a systematic error which shows up as a non-zero 
AH = 0(At 2 ). The detailed- balance is not endangered by this error, because a 
Metropolis acceptance step is performed. For fixed trajectory length, N can then 
be tuned to optimize the overall efficiency. 

To generalize the method to constrained variables, one needs to appropriately define 
the conjugate momenta and the equations of motion in order to preserve the constraint 
and, most importantly, not to spoil the time reversibility. Each spin variable, 0, lives on 
the surface of a two-sphere, and correspondingly one could imagine an algorithm with two 
independent conjugate momenta, maybe living in the perpendicular plane (0 • 7r = 0). 
However, changing the constraint from the field to the momenta is not very appealing 
(and, from the practical side, one would need to worry about two constraints in the numer- 
ical integration). A different approach, the use of spherical coordinates, has the drawback 
of a non-planar integration measure. Our very simple algorithm avoids constraints and 
non-planar measures, by introducing three conjugate momenta per spin. 

We shall start from an analogy with the dynamics of a particle living in the sphere, a 
potential (V) acting on it. The Hamiltonian is 

^sphere = — + V((f)). (40) 

Here L is the angular momentum, 0x0. The equations of motion can now be obtained 
from the Poisson Bracket |37|] with the hamiltonian (|40|): 

SV 

= Lx0,L = -0x— . (41) 

SV j„ r_ / SV 8V 5V ' 



In this expression ^ stands for ^ 
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This formalism is still inconvenient for us, because the constraint cf) L = complicates 
the generation of random momenta according to a Gaussian distribution. However, the 
following simple facts can be straightforwardly established from the equations (|4l|): 

I. Both cj) 2 and cf> ■ L are conserved through the time evolution. If the initial condition 
verifies the constraints 4> 2 = 1 , 4> ■ L = 0, this will not be spoiled by the dynamics. 

II. The dynamics is time-reversible. 

III. Although the Li cannot be all canonical variables ||37|| , the "Liouville" measure, 
D(f) DL{= dcfii dcp2 d(j)3 dL\ dL2 dL^), is left invariant by the time-evolution. 

IV. The Hamiltonian is a constant of the motion. 

Now let us forget about the constraint (j> L = 0, i.e. we introduce a new field P which can 
have a "radial component" (it is no longer an angular momentum), but we keep the Eqs. 
of motion (|4~l|). Obviously, statements I-IV will still hold. Whether a symplectic structure 
is hidden under this new dynamical system is unclear, but also irrelevant (properties II 
and III are the essential ones for HMC to be a correct algorithm p6[). 

So, we introduce three momenta per spin, P = (Pi, P2, P3), and write down the 
Hamiltonian 2 

X 

Equations of motion respecting properties I-IV are easily generalized: 

S S 

4>(x,t) = P (x,t) X 4>{x,t) , P (x,t) = -4>(x,t) X TT • ( 4 3) 

( "r (x,t) 

As expected, the evolution equations for the S 2 fields take the form of (infinitesimal) 
rotations, while the conjugate momenta can be considered as living in the Lie Algebra of 
SO (3). The discretized leap-frog form of these equations is therefore naturally formulated 
in terms of finite SO(3) rotations, 

^(nAr + Ar) = exp[ArP x {(n + -)Ar) • J] ^(nAr) , (44) 



P x ((n+l)Ar) = P x ((n - l)Ar) - ( ,, nAr) x ? Sb At , (lo) 



where J are the 3x3 generators of SO (3), satisfying 

(exp[#n • J])ij = 5ij cos 8 + riirij (1 — cos#) — e^kfik sin 9 (46) 

for unit vectors n. Again, the length constraint on the <f> fields is preserved by construction. 

This final result is reminiscent of the elegant solution for models with variables be- 
longing to a Lie group and conjugate momenta in the group algebra (or vice versa) p5| . 

In our case, Sfiquenched = ~k^2 nfl 4> n ' ^n+jii so the HMC algorithm can now be 
implemented in a straightforward manner. To test the algorithm, we have simulated 
the 0(3) model on an 8 3 lattice at k = 0.693 k c @ with our HMC algorithm and 



with Wolff's single-cluster embedding algorithm [39]. Let us first define the measured 
observables, and then compare them. 

In this work we have only measured bosonic observables, as our sole objective was the 
numerical determination of the phase diagram. We have constructed our observables in 
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Table 1: Values for several observables in the quenched model (|3]) on an 8 3 lattice at 
k = 0.693 ~ k c , obtained with our implementation of HMC and with Wolff's single 
cluster algorithm p9| . 



Algorithm 


(E) 


d k (E) 


x/v 


£ B 


HMC 


0.3505(5) 


1.51(2) 


0.1426(9) 


4.47(2) 0.800(6) 


Wolff 


0.35061(13) 


1.501(10) 


0.1432(2) 


4.486(9) 0.8031(18) 



terms of the Fourier transform of the spin field: 



£=o~o V , ( 4 9) 



as 

where V = L 3 is the lattice volume. 

We define the non-connected finite-volume susceptibilities as 

x = V(m 2 (0,0,0)>, Xs = V(m 2 (TT,7r,TT)). (48) 

The subscript 's' on Xs stands for 'staggered', and this term is used to label quantities 
which are taken with a weight —1 for the odd lattice sites, corresponding to momentum 
(tt, it, 7r). Notice that x/V is a pseudo order parameter, which should be of order one in a 
ferromagnetically broken phase, and of order 1/V in a paramagnetic or antiferromagnetic 
phase (and similarly for x s /V). 

Another quantity of interest is the Binder cumulant 

5 3 ((m 2 (0,0,0))' 
2 ~ 2 (m 2 (0,0,0)) : 

with an analogous definition for the staggered variant B s . 

One expects S = 1 in the FM phase, where x/V is non- vanishing in the thermodynamic 
limit, while it should be of order 1/V in the PM phase, far from the phase transition. 

For the correlation length, we use a definition which is easy to measure and gives 
accurate results: 

X/F-l V /2 
v 4sin 2 

where F is the squared Fourier transform at minimal non-zero momentum, 

F = — ^|m(2vr/L,0,0)| 2 ^> + permutations) . (51) 

Again, the generalization to staggered quantities is straightforward. Another kind of 
observable, needed for the standard extrapolation method [Q| , is the normalized nearest- 
neighbour energy 

E = ^E^'^+a> = Yk^ z - (52) 
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We also measure its fluctuation, given by 

3V (<£ 2 } - (E) 2 ) = JL (E) . (53) 

In Table |l] we compare the values obtained for these observables, using our HMC 
algorithm and the single-cluster algorithm. We find excellent agreement. Of course the 
efficiency of our implementation of HMC is not competitive with a cluster method in the 
0(3) non-linear a- model. But it could be useful in other models where cluster methods 
are not effective in reducing the dynamical critical exponent z (for instance, when some 
kind of frustration is present ^0|), while HMC is expected to yield z = 1 for any bosonic 
model. 



4.2 The HMC algorithm for the full theory 



The only restriction imposed on HMC is that the fermion bilinear in the action should be 
given in terms of a positive definite matrix. This will the case if we consider two identical 
fermion families (Nf = 2) as is usually done in lattice gauge theories. After integrating 
them out we obtain (detM) 2 = det(M'M), where M is the fermion matrix for a single 
fermion family. As we are mainly interested in the strong spin-fermion coupling region, it 
makes sense to perform the following manipulation: 

det M = det (Y + K) = y AV det (1 + Y~ l K) (54) 

(cf. Eqs. (30,31])). The constant factor y 4V can be dropped, and we define M = l + Y^K . 

Next, one re-exponentiates the (inverse) fermion matrix by introducing the so-called 
pseudo-fermions z x , which are complex four-component c-number fields. The partition 
function is then 

(55) 



Z = JD(j)DzDz exp (-£* B - ^(M+M)- 1 ,? 



For further details we refer to Ref. [36]. 
Now the HMC Hamiltonian becomes 



H 



1 



x+li + z ] (M f M 



and the time reversible, constraint and energy preserving equations of motion are 



4>(x,t) 



(x,t) 



(x,t) 




(56) 



(57) 



M + h.c. 



-1 



For the inversion of the fermionic matrix, we have employed the conjugate gradient 
algorithm. To formulate the stopping criterium, let us define h = (M' Af) z, h n being 
the n th trial solution. We continued the conjugate gradient iteration until 

2 

< R. (58) 



\(M*M)hn - z\ 



\hr. 
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In the simulation, we need the inverse matrix both for the leap-frog and for the 
Metropolis accept-reject step. It is clear that R does not need to be the same in both 
cases. For the Metropolis step, lack of accuracy in the inversion will bias the simula- 
tion. To control this, we have checked that the Creutz parameter (exp(-AH)) equals 1 
within errors. In some regions of parameter space R values as small as 10 -25 were needed. 
The essential requirement on the leap-frog is full reversibility in the numerical integration 
of the equations of motion (up to the numerical precision reachable with 64-bit floating 



point arithmetic). As first noticed in ref. 41], this has no relation with R if the seed 



for the conjugate- gradient algorithm is chosen to depend on the actual configuration only 
(ho = z, for instance). However, if R is too large, the numerical integration will produce 
large changes in the Hamiltonian, and the Metropolis acceptance will be poor. We have 
found that R = 10~ 7 during the leap-frog steps allows for a 50% acceptance. 

In a first implementation of a new MC algorithm, some consistency checks are ex- 
tremely useful. In addition, there are three parameters to be adjusted for optimal perfor- 
mance, N,t and R. We have carried out the following tests: 

1. We have explicitly checked reversibility of the leap-frog algorithm. 

2. We have checked that (exp(— AH)) = 1 within errors. 

3. The gaussian expectation values, (M^M) 1 = 4 and (P 2 ) = 3 have been 
checked. 

4. We have checked that AH oc (At) 2 in the leap-frog integration, for constant trajec- 
tory length NAt. 

In addition, we compared simulation results for the full theory at (k,y), with the 
output of a quenched simulation at the corresponding effective coupling value obtained in 
a large- y expansion, 

y \y 

(cf. Eq. (pi])). In table (|2|), we give the mean value of several operators as obtained 
on a 4 3 lattice at k = 0.693, y = 10.0 and in the quenched theory. The agreement is 
excellent. Notice that even if the shift in the effective coupling is only 3%, the effects of the 
dynamical fermions can be clearly measured as the observables change quite significantly 
at the critical point k c = 0.693. 



fceff = k + 72 + ° ( 4 ) ( 59 ) 



4.3 Phase Diagram 

The phase diagram in fig. [I] was obtained on an 8 3 lattice. As there is no true phase 
transition on a finite lattice, a criterium is needed to locate the phase boundaries. We 
looked for the point where the relevant Binder cumulant equals the value B = 0.8 it has 
at (k = ±0.693 ~ k c , y = 0). Since B = 1 deep in the broken phase and B oc 1/L 3 in the 
symmetric one, this provides a clean quantitative criterium which yields a point definitely 
inside the critical region. The width of the critical region decreases as L^ 1 ^ , therefore the 
systematic error in the critical coupling will be at most of order 10 _1 . However, since the 
Binder parameter is a universal quantity, which should stay constant along much of the 
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critical lines, the error rather goes as 

jr-( w +i/„) ^ O ( 10 -2)). Thus, this systematic error 
is under control in the full theory as well. We used the Ferrenberg-Swendsen extrapolation 
method to determine the precise location of the points where B = 0.8. 

The total simulation time was 16 days of the 32 Pentium Pro processor parallel com- 
puter RTNN based in Zaragoza. To allow for a correct thermalization, we discarded 100 
integrated autocorrelation times of the relevant susceptibility. This may look utterly con- 
servative, and the MC history indeed seems to stabilize long before that. However, not 
much is known about the exponential autocorrelation time of fermionic algorithms and 
one should be cautious. 

As Eq. shows, both at y = oo and at y = we recover the non-linear cr-model 
with its well known paramagnetic, ferromagnetic and antiferromagnetic phases. At finite 
y, we expect these phases to extend into the (k,y) plane. In fact one can quite precisely 
anticipate the critical coupling from the strong coupling formula (59) and the quenched 
critical points k^ -00 ^ = ±0.693. Using the Ferrenberg-Swendsen extrapolation procedure, 
the phase transition lines can be determined down to y ~ 2.0. In fig. ^ the variation of 
the Binder cumulant and the susceptibility around the two critical couplings is shown for 
.</ 2.1). 

In the small-?/ region, the effective action up to 0(y 2 ) does not only renormalize k, but 
also introduces additional couplings, due to the non-locality of the matrix K~ l occurring 
in the weak-coupling expansion. Therefore, we do not have an estimate for k as reliable 
as in the large- y region (|5^), but we can nevertheless obtain an estimate for k c (y) from 
the MF approximation. We have simulated at several values of the coupling k, for fixed y, 
until the corresponding Binder parameter crossed its critical value. A more accurate result 
for the critical point was later on obtained with the Ferrenberg-Swendsen extrapolation. 
In fig. ^, we have plotted the relevant Binder parameter and susceptibility for k values 
near the two critical couplings with y = 0.5. 

In fig. ^| we show the variation of both order parameters and Binder cumulants when 
crossing the FM(S)-EM transition line at y = 1.15. We find a strong change in the stag- 
gered quantities, while the non-staggered ones show a smoother evolution. However, the 
non-staggered order parameter is much smaller than its staggered counterpart. This may 
indicate that, although the non-staggered sector is non-critical {B ~ 1), it will eventually 
undergo a phase transition at lower k. A similar behaviour is found when traversing the 
AFM(W)-EM line at k = —1.6 (see fig. |5|), but now the non-staggered quantities show a 



Table 2: Comparison of observables in the full theory @ at (k = 0.693, y = 10.0) 
and in the quenched model both at the corresponding value of k eS and at k c = 
0.693. We have 140, 000 unquenched trajectories (i\/=10, Ar=0.3) on a 4 3 lattice. 
The Metropolis acceptance rate was 65-70%, with an autocorrelation time of 3-4 
trajectories. 



Couplings 




(E) 


d k (E) 


x/v 




fc=0.693 , y= 
k=0.7U , y= 
A;=0.693 , y= 


=10.0 

=0 

=0 


0.4164(6) 

0.41584(14) 

0.3928(3) 


1.134(6) 
1.130(4) 
1.174(4) 


0.3111(7) 
0.3108(2) 
0.2836(4) 


2.378(4) 

2.3779(18) 

2.214(2) 
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Figure 2: Binder cumulant ( |49| ) and non-connected susceptibility (^) as a function 
of k, around the two critical points at y — 2.0. For each critical point, only one 
simulation has been carried out. The other points are obtained with the Ferrenberg- 
Swendsen extrapolation method. 



more pronounced signal. The detailed study of these transition lines (order of the phase 
transitions, critical exponents, etc.) requires a finite-size scaling analysis, which is left for 
future work. This study will be much easier if the transition line is crossed varying k, as 
we lack an analogue of the Ferrenberg-Swendsen extrapolation method for y. 

5 Quasiparticle excitations at the MF level 

In this section we explore the relevant excitations involving fermions, with emphasis on 
the strong-coupling region of our model. This will then enable us to discuss electronic 
properties. 

The small-?/ regime has been studied in relation with the mechanism by which leptons 
and quarks acquire their mass through symmetry breaking in the Electroweak sector of the 
Standard Model. Due to the weak coupling, one has essentially Fermi liquid behaviour, 
and there are no surprises. This situation will change dramatically when we consider the 
strong-coupling region, though. 
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Figure 3: Binder cumulant ( |49| ) and non-connected susceptibility (^) as a function 
of k, around the two critical points at y = 0.5. The data points are from different 
simulations. 

5.1 Fermionic excitations in the FM(S) and PMS phases 

At very large y, it is natural to attempt a large-y expansion. This can be achieved after 
carrying out the following change of variables: 

& = (60) 
= (0-t)V>. (61) 

Because of the constraint 4> 2 = 1 and the identity (cf) ■ r) 2 = (f> 2 l, this transformation has 
unit Jacobian and its inverse satisfies 

$ = (0 • r) V' ■ (62) 
In terms of the new variables (dropping the primes) the action takes the form 

S = -k ^2 <t> x ■ 4>x+? + Yl $ x ( Kx v (&V ' T ) + y6x v) ( 63 ) 

x,n x,y 

where the fermion kinetic term is the usual lattice kinetic Dirac operator, defined in Eq. 
(|l2|). After a further rescaling of the ip fields, the coupling y can be moved to the kinetic 
term, where it appears as 1/y. 

Note that this change of variables (^J6^) was implicitly present in the MF calculations 
of the phase diagram in the strong-coupling region as well (cf. Eqs. ( |30| , [3l| )). 
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The fermion propagator {ipx4>y) 1S given by the expectation value of the inverse fermion 
matrix, which in a large-y expansion becomes 

{il> x $ v ) = (M- 1 ) = ( - (l - -K{<t> ■ t) + \k(4> • t)K{4> ■ t) - • • •) ) . (64) 

\ y v y y j X y I 

This can be viewed as a sum over paths of increasing length connecting x and y (recall 
that if is a nearest-neighbour matrix). 

In the FM(S) phase, there is a non-zero magnetization v = \ {4>)\- Expectation values 
of products of (f> fields on different sites are replaced by the appropriate powers of v. 
Corrections to this approximation as well as contributions from paths visiting a given site 
more than once are of higher order in 1/d and are ignored at the MF level. The series 
(|64|) can thus be resummed and one finds a propagator 

1/v 

W)fm{S) = K + y/v ( 65 ) 

which is that of a fermion with a mass y/v. Note that, since v < 1, this is a huge mass if 
y is large. The propagator for the original fermion, before the change of variables ( |60|j6l| ), 
corresponds to the same physical particle; the only difference is in the wave function 
renormalization. 

In the PM(S) phase, v = 0, so at the MF level the fermion would be infinitely massive, 
or in other words, non-propagating. Beyond this naive MF level, however, a large but 
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Figure 5: Binder cumulants and susceptibilities when crossing the AFM(W)-EM 
transition line at k — —1.6. 



finite mass will be found. This is due to the next-to-leading contributions to the series 
(p4[). The dominant terms are now those involving the expectation value for the nearest- 
neighbour energy z 2 = (cf> x ■ 4> x+ ^}, which is of order l/2d and therefore absent at the MF 
level. The resummation of contributions in (^34|) now leads to a fermion propagator with 
a mass y/z, which is even larger than the mass of the fermion in the FM(S) phase. 

The conclusion of this analysis, which is similar to that in (chiral) Yukawa models 
in the Electroweak theory |44|] , is that the elementary fermion excitations in the large-y 
region are very heavy (hence essentially non-propagating), and therefore play no role in 
the spectrum of light excitations. This holds even stronger in the PMS phase than in the 
FM(S) phase. 



5.2 Fermionic excitations in the AFM(S) phase 

Here our point of departure is again the form of the action (|63|), which is tailored for 
studying the large-y behaviour. In the AFM(S) phase, we have a staggered expectation 
value for the cf) field at the MF level, which can be taken in the 3-direction, 



4> x = e x v\ (66) 
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(with e x = (_i)*i+*a+*s). Hence 

(0 X • T) 



(67) 



where ipi\ i = 1,2 labels the two flavours in ?/V So after the change of variables ([30.61) 
the kinetic operator in (^) is still diagonal in flavour. The only effect of the new variables 
is to change the lattice Dirac operator from (^) to 



v e y t 3 K, 



Due to this diagonal structure in flavour space, we can concentrate on one flavour, say 
^W; the other flavour is obtained by taking — v instead of v. In Fourier space, the kinetic 
term for ip^ 1 ' is given by 

- i v s/np 5 Ptq±7r , (68) 



where 



s/np 



(69) 
(70) 



So we obtain for the inverse of the MF propagator in the AFM(S) phase, 

M p , q = -i v s/np <5 Pj(?±7r + y S Pjq , (71) 
or, in matrix notation for the subspace of the modes coupled in Eq. ([71]), p and p±(7r, vr, 7r), 



Af. 



p,pi(7T,7r,7r) 



y — i vsfoip 
i vsfoip y 



(72) 



To find the quasiparticle excitations in the AFM(S) phase we diagonalize the fermionic 
part of the action (|72j) . One obtains 



S = ip{p) (y — v s/np) ip(p) 



(73) 



where 



or, in position space, 



V2 



ip W (p) + i^ w {p + ir) 



l 

75 



+ie x ^ 



The momentum space propagator corresponding to (|73j ) is thus 

1 y + v s/n p 



S(p) 



y — v s/np y 2 — v 2 J2\ sm2 Pa 



(74) 
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Since we are working in imaginary time, one would expect quasiparticle poles in S(p) to 
appear at negative values of p 2 . The unusual relative minus sign in the denominator ( |74| ) 
therefore does not seem to allow for a quasiparticle interpretation, at first sight. 



However, (74) suggests the possibility of light excitations with a relativistic dispersion 
relation around spatial momenta (±^,±^). To see this, consider the denominator in Eq. 
(f74|) for small = p^ ± ir/2: 

y 2 -v 2 Y, PA = ~ v 2 d) + v 2 k l + ^(k 4 ), (75) 
A A 

where d = 3 is the space-time dimension. As long as we are at large enough y, such that 
y 2 > dv 2 (recall v 2 < 1), this dispersion relation corresponds to a relativistic excitation 
with m 2 = (y 2 — dv 2 )/v 2 , in this naive MF calculation. Several comments are in order: 



1. For v=0, we recover the MF result for the PMS phase: the kinetic term in (|75| ) is 
suppressed. 

2. At the MF level, only for (y 2 — dv 2 ) small enough compared to v 2 these fermionic 
excitations, (r • <fi) ip, can propagate easily. Since v 2 < 1, this can only happen for 
y 2 not too large. Whether or not this situation will arise depends on the precise phe- 
nomenological relation between x and (y, k), which is likely to be dopant-dependent 



(see section 2.2). 



3. These would-be excitations are characteristic of the AFM(S) phase. Let us recall 
that in the PMS phase no light fermionic excitations have been identified at the 
MF level. In any case, in the SC phase of our model we are not in an AFM 
background. Therefore no momentum-space rotation will be needed to interpret the 
fermionic excitations and no hole-pockets around (±7r/2, ±7t/2) are expected in the 



SC phase g8|. 



4. We need to fix the scale in some way, as done for instance in 0], to estimate by 
means of a MC calculation the masses of all possible excitations of our model. 



5.3 Light bound states in the PMS phase 

We have seen above that the fermionic excitations in the PMS phase are very heavy. We 
will now show that there are bound states of elementary fermions, however, which are 
light. This is done by means of a MF calculation of the double-chain type These 
pairs will be the electron or hole pairs responsible for the superconductivity, as we will 
discuss in Sect. [6|. 

Consider the propagator for a fermion pair ^ x ip X : 

(^4%) = ( M *,kyXk M *li;yJ ~ ( M *,li;y,X,k M ^kr,yJ " i 7(i! 

Here M _1 is the single-fermion propagator, a,(3,\,p are Dirac indices, and i,j,k,l are 
flavour indices. Thus, this propagator is really a 16 x 16 matrix. For the moment we keep 
all these indices as they are; later on we will discuss how pairs of them decompose into 
quantum numbers for the composite state. 
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Figure 6: A typical double-chain diagram, connecting sites x and y. The chains 
are parallel in position space. 



Let us concentrate on the first (M 
of M _1 as before, we find the series 



M~ 



term in Eq. (|7q). Using the 1/y expansion 
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y V v 



N 
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x,a,i; y,p,l I 



(77) 

where we have written 4> as a shorthand for (</> • r). It is clear that only terms with 
N + N' even survive in a paramagnetic phase, due to the cf> — ► —<f> symmetry, thus a 
factor (-1)^+^' has been dropped. Since the matrix K connects nearest-neighbour sites 
only, each term in this series can be seen to represent a product of two paths (chains) of 
lengths N and N' respectively, connecting site x with site y (so, if the "distance" between 
x and y is even(odd), both N and N' will be even(odd)). 

We will attempt to sum the complete series, to leading order in 1/d, where d = 1+2 = 3 
is the euclidean space-time dimension. For this, we need the spin-spin propagator, which 
in this approximation is extremely short ranged 



1 



■ab 



(78) 



Expectation values of the type (4> x -4> x+ jx) are of order 1/d, and others are suppressed even 
stronger. Thus, assuming (|78|), we observe that any term in the series which contains <f> x 
for a given site x only once or an odd number of times will vanish due to ((f)) = 0. When 
the site is visited twice, it follows from 2 = 1 that the contribution from the fields is 
proportional to \& ab ■ Thus each site along the chains connecting x and y must be visited 
an even number of times. One class of diagrams fulfilling this requirement consists of the 
so-called 'double-chain' diagrams, where the propagation of both fermions between x and 
y follows the same path in position space (see figure As was convincingly argued in 
Ref. [43], this class saturates the dominant diagrams in the 1/d expansion. Indeed, one 
can easily check by concrete examples, how deviations from double-chain behaviour induce 
additional powers of 1/d. We shall also assume that these double chains are self-avoiding 
(this is allowed at first order in 1/d). 
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Figure 7: A matrix-product term contributing to the flavour structure. 



Our task is thus to sum up all double chain diagrams connecting x and y. Let us first 
consider the flavour structure. Using (78) one finds that 



{(<f>x- T )jk(<f> x -T)ii) = \^2r°- k T°i = -(-S jk Sii +2SjiS ik ) 



(79) 



From this and from the ultra- local correlations we are considering (cf. Eq. (|78j)), it follows 
that the product of 2{N + 1) factors of (0 • t) along a double chain of length N visiting 
the points x = xq,xi, . . . ,y = xjy (cf. Eq. (|77|)) is 
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(80) 



x,i;y,l' 



To calculate P and Q, it is convenient to represent the general term contributing to the 
above matrix product as in figure |?]. A graph contributing to 5jk5u will have an even 
number of crossings, while diagrams contributing to SjiSik jump an odd number of times. 
Each crossing contributes a factor |, while non-crossings yield factors -§ (cf. Eq. (j7|)). 
Now, P and Q can be easily obtained using binomial summation: 
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■^(Sjkdii + SjiSik) + (-1) -{5jk5u - SjiSik) , (81) 



where we have separated in a term symmetric under (ji) *-* (ij) and an antisymmetric 
one (this will be needed for separating the contribution to different quantum numbers). It 
is remarkable that the flavour contribution only depends on the double-chain length, but 
not on its shape. This allows for a total factorization between flavour and Dirac indices. 
Next, consider the Dirac structure. One gets products of matrices 



along the double chain, where 



(82) 



(83) 



One readily finds that 



Axy — ^ -^xy-^xy ~ {^y,x+jl + 3y,x-fi) • 



(84) 
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Thus, we need to calculate 
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(85) 



x,/3,a;y,X,p 



where the sum is extended to all the lattice paths (denoted by {fJ. n }) of- length N connecting 
x and y. Now, we can extend the sum to all length- A lattice paths starting at x, because 
paths not connecting x to y will contribute a zero x y entry. This can be also understood by 
realizing that once the chain has arrived at Xj, there are 2d possible directions to continue 
the chain. These are added up by summing Eq. (]8^) over fi. At the next site, we do the 
same for the next step along the chain. The contributions of all double chains are therefore 
added up when we take the product of these /i-sums along the chain. Corrections due to 
backtracking (2d — > 2d — 1) are down by \/d. 
So we need to calculate powers of the matrix 



(86) 



One way to do that is to write it out explicitly as a 4 x 4 matrix in the space spanned by 
the vectors (/?, A) =(1,1), (2,2), (1,2) and (2,1), in that order. One finds that it equals 
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It can be diagonalized in this 4x4 space. The eigenvalues, up to the factor 1/4, are found 
to be 



A" 
A 4 



A - 2A" 
-A, 



(^ = 1,2,3), 



where 
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□ + 2d, 



(89) 



(90) 



and □ is the lattice discretization of the d'Alembertian ^Z^d^d^. The A power (see 
flS5"D) of the matrix ( |S6"| ) is now easy to calculate. 

In order to collect the factors and sum up the contributions, let us go back to Eq. 
(f76|). We see that we need to antisymmetrize each term in (Af with respect to the 

simultaneous interchange of a,i with (3,j. This gives a sum of two terms, one symmetric 
in a <-> (5 and antisymmetric in i j (corresponding to a composite state which is a 
Dirac vector and a flavour singlet), and one vice versa (singlet in Dirac space, vector 
in flavour space). Note that Eq. ( |sl| ) has already been written as a sum of symmetric 
and antisymmetric terms. The symmetric and antisymmetric parts of the Dirac structure 
correspond to Eqs. (|88|) and (|89|), respectively. 

Collecting the various factors, we can carry out the geometric sum over in Eq. ( |7"6| ) 
and we find the following propagators for the composite states: 



28 



• a Dirac vector - flavour singlet with propagator 

-□ + 2A» - 4y 2 - 2d 
where a, b are the Dirac vector indices 

• a Dirac singlet - flavour vector with propagator 

-85 u 
-□ - 12y 2 - 2d 

where /, J are the flavour vector indices. 

These have the form of massive bosonic propagators, up to the following caveat (of 
course, higher order corrections in 1/d may induce shifts in the precise location of the 
poles, as well as their residues). 

The propagators in ([)!]) contain the matrix 2A^ in the denominator. However, this 
term must be ignored since it is sub-dominant in 1/d, compared with the (lattice) d'Alembertian 
□ . 

The numerator of the propagator (^) carries a delta function only, instead of the usual 
tensor structure — d^d u /m 2 . This is also an artifact of the 1/d approximation. 

Notice also that the terms which would play the role of a mass squared in the denom- 
inators have an apparently wrong sign. However, it is easy to check that the composite 
field e x if) x ip x (where e x = (— l)^ 2 ^ as usual) does lead to a massive Dirac singlet - flavour 
vector propagator with mass squared m 2 Q ^ = 12y 2 — 2d= 12y 2 — 6. Similarly, one obtains 

a massive Dirac vector - flavour singlet with a mass squared ^ = 4y 2 — 2d = Ay 2 — 6. 
We thus conclude that the right interpolating field is e x ip x ip x |4*| (one could equally well 
consider it as an excitation centered around spatial momentum (vr,7r), though). 

The conclusion is that we find massive bound states of fermions in the PMS phase. 
They are bound by the strong interactions with the spin waves. These composites are 
lighter than the elementary fermions in this phase, when y moves away from the value do. 
This means that they will be the dominant light excitations, and we will argue that they 
condense and lead to superconductivity in our model. 

6 The x—T phase diagram of the cuprates 

In this section we shall consider the application of the above results to the cuprates. 
Strictly speaking, our calculations only apply to the zero-temperature case. However, 
some conjectures can be made for the case of non-zero temperature, by means of very 
general thermodynamic considerations. 

6.1 Physics at zero temperature 

We will discuss how the copper oxide materials are described by our model, by following 
the 'trajectory' which they map out in the y, k phase diagram of our model as the doping 
fraction x is increased from zero all the way to the overdoped regime. Along the way, 



(91) 



(92) 
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various phases will be traversed, and the relevant light excitations present will be given 
their physical interpretation. 

As we discussed earlier on, it is by now well established that the undoped material can 
be described by a point (k < k c ,y = oo) in the AFM(S) phase of the 0(3) model to which 
our model reduces for y = oo. The actual value of k is related to the spin stiffness and 
velocity; in the large-5 approximation k oc S 0. As the doping fraction x is increased, the 
carrier mobility increases which we have argued in section 2.2 to correspond to decreasing 
y in our model. 

We have schematically indicated this 'evolution' of the cuprates with increasing doping 
by the line of arrows in Fig. |l[ For illustrative purpose, let us assume (see our comment 
in section |2.2| ) a relation of the type x ~ C 2 /(C 2 + y 2 ) for some constant C (in the case 
C = 1 the horizontal axis in Fig. [I] would then correspond to 1 — x), but this is immaterial 
for the qualitative picture. 

Note that we move in the direction of the PMS-AFM(S) transition line, so the AFM 
order will decrease. This is consistent with the experimentally observed reduction of AFM 
order upon doping. 

At some point along our trajectory, still within the AFM(S) phase, the possible exci- 
tations at the (±7r/2, ±7r/2) hole pockets (cf. Sect. f>.2\) would become light and start to 
dominate. 

When the doping is increased even more, we move into the PMS phase. As long as we 
remain close to the PMS-AFM(S) transition in this phase, short-range AFM correlations 
will still be present, although they are predicted to decrease as one moves deeper into this 
phase. The crucial point is that, as demonstrated in the MF calculation of Sect. |5.3| , the 
only light excitations left now are the "PMS pairs", bosonic fermion bound states. Let 
us remark that a rich spectrum of excitations with different quantum numbers (not all of 
them identified in our MF calculation) should be expected. In our calculation, the flavour 
singlet (i.e. the physical spin-singlet) has turned out to be lighter than the triplet, which 
is encouraging. To investigate this fundamental point in more detail, a MC calculation 
has to be done. 

At temperature T = these pairs will be Bose-Einstein (BE) condensed, as any other 
bosonic state would. This leads to superfluid behaviour for the pairs, and to superconduc- 
tivity once electromagnetism is coupled into the model. Since a finite number 0(x Cl /2) 
of them becomes available at the same time, at the point where we enter the PMS phase 
(corresponding to some critical doping x Cl ), one expects a finite (i.e., not infinitesimally 
small) critical temperature T c = Tbe here. 

Following the arrow line towards even smaller y, we leave the superconducting PMS 
phase at x C2 . Thus we understand why superconductivity only happens for a range of 
doping. Again, T c is expected to remain finite up to x C2 . This may explain why T c is 
experimentally observed to jump steeply at x Cl2 ; however, we will discuss the possibility 
of even more unusual behaviour when we discuss the behaviour below. 

Our model predicts a ferromagnetically ordered phase at T = in this intermediate-?/ 
regime. The system traverses a region where one needs to go over from a description of 
the charge carriers in terms of the variables (r • 4>)ip at large y, to a description in terms 
of the fields tp appropriate for the small-y Fermi-liquid regime. (There may be a kind of 
cross-over in this region.) In either description, propagation of the carriers is hampered 
by the effective mass induced by the FM spin waves. 
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At this point one may object that no FM behaviour has been observed experimentally 
in the overdoped regime. We note, however, that at some non-zero temperature (cf. Sect. 
p^ ), this phase will be replaced with a thermally disordered phase. Since the overdoped 
regime has not been explored in as much detail as the underdoped and superconducting 
regimes, in particular at very low temperatures, this may explain why no FM behaviour 
has been detected here. 

For very large doping we reach the weak-coupling region with Fermi-liquid behaviour. 
We will end up in the PMW, AFM(W) phase or possibly remain in the FM(W) phase of 
our model. This will depend on the actual trajectory of the system with doping (see our 
comments in section |2.2| , and recall that for large doping fraction the effective k parameter 
might also be influenced by it). 

All these considerations illustrate how such widely varying behaviour in the cuprates, 
controlled by the doping fraction x, is reproduced by varying just one parameter in our 
simple two-parameter model. A point needing further investigation is the following. At 
strong coupling, one would expect a shift from commensurate to incommensurate AFM 
ordering, with increasing doping |4q| . After integration out the fermion field in our model, 
frustrating couplings (of order 1/y 4 in the strong coupling expansion) are generated, which 
presumably lead to such a phenomenon. Our MC simulation, was done on a too small 
lattice to be able to resolve such an effect, and our MF is not reliable in this intermediate-?/ 
region f33|| . 

An interesting prediction of our model is that superconductivity is unlikely to occur 
in materials with spin S > 1/2. The reason is that the undoped model would correspond 
to a point k <C k c in the 0(3) model (\k\ oc S, see ref. [[7|), making it unlikely that 
the "evolution trajectory" pass through the PMS phase. In more physical terms, the 
larger \k\ is in our model, the stiffer is the magnetic ordering and the less effective are 
the spin fluctuations needed to bind the fermions. For instance, upon doping the layered 
compounds Lai-zSri+^MnO^ which have localized S = 3/2 spins (implying k ~ — 2 
in the 0(3) model) a disappearance of the antiferromagnetic phase and the subsequent 
emergence of an exotic magnetic phase (maybe a spin-glass phase) is observed |46[| , but no 
superconductivity appears. Essentially, the same thing happens for insulating nickelates 
L^-zSr^NiO^ for which the localized spin is S = 1 [47|. 



6.2 Physics at non-zero temperature 

What happens when the temperature is increased from zero? So far, we do not have 
any MF or MC results at non-zero temperature. However, by the following fairly general 
arguments we are led to conjecture Fig. |8| as a schematic sketch of the most general x-T 
phase diagram suggested by our model. 

First, let us discuss what happens with the superconducting phase when the tempera- 
ture is increased. Consider a doping fraction x corresponding to a y value such that we are 
well inside the PMS phase, with superconductivity at T = 0. At a certain temperature, 
Tbe, the BE condensation will be undone, so the superconductivity disappears. We thus 
identify Tbe as the critical temperature T c . Above T c the fermions are still bound together 
in (uncondensed) light bosonic PMS pairs. One therefore expects "normal", charge-2 con- 
duction due to these pairs. This phase is sometimes called a "spin-gap" phase. As a matter 
of fact, the existence of such a phase is shared by any model in which superconductivity is 
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Figure 8: Sketch of the predicted phase structure in the x-T plane. The dotted 
vertical line at xs-w indicates a possible crossover. 

triggered by the Bose-Einstein condensation of previously formed fermion-fermion bound 
states. Experimental evidence for such a mechanism is now available ||. 

Above some spin-gap temperature, Tsg, the PMS pairs will be broken apart by ther- 
mal fluctuations. However, notice that the pair-breaking temperature for the first, isolated 
PMS-pair is lower than the temperatures needed to dissociate many of them, since the 
constituent fermions have to fill up states of increasing energy in the Fermi sea. Thus one 
expects Tsg to be characteristic of a cross-over rather than a well-defined phase transi- 
tion. In this high-temperature phase, the only possible carriers would be fermions, but 
as we have already discussed, their kinetic term vanishes in the mean field approxima- 
tion. (This perfectly insulating behaviour will receive corrections beyond the mean-field 
approximation, though.) 

In the qualitative behaviour shown in Fig. |8|, we have supposed that Tsg is always 
larger than Tbe- Depending on the variation of the mass and the binding energy of the 
bound state with doping, it could happen that the Tsg and Tbe critical lines join at 
some point. The actual situation regarding the full curves Tsg and Tbe would have to be 
investigated numerically in our model. 

Further qualitative statements can be derived from the general thermodynamic ar- 
gument that raising the temperature has the effect of increasing the magnetic disorder. 
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In the phase diagram of our model, Fig. |l], this means that the two PM phases expand 
in all directions, at the expense of the various magnetically ordered phases. One of the 
consequences could be a rash disappearance of the FM phase, which might explain why it 



has not be observed, as we briefly discussed in Sect. 6.1 



Another consequence is the intriguing possibility of re-entrant superconductivity for 
certain values of x, with superconductivity not starting at K, but restricted to a temper- 
ature interval < Tm <T <T c . To see this, consider x slightly smaller than x Cl , corre- 
sponding to a point (y, k) in figure [l] inside the AFM(S) phase, very close to the transition 
to PMS. Let us assume the disorder due to the temperature has similar properties to the 
disorder coming from the dynamics (i.e. at zero temperature). Then the PMS phase will 
be enlarged and absorb the point (y, k), implying superconductivity provided that we are 
still below Tbe- Similar re-entrant behaviour is expected for x slightly larger than x C2 . In 
order to make quantitative predictions, an estimate of the mass and the binding energy 
of the PMS-pairs would be required. Some experimental evidence for reentrance has been 
found ten years ago jl8| in YBaCuO, but it seems to have gone largely unnoticed. Further 
experiments are required to resolve this issue. 

In fact, our model does not exclude the extreme case in which the superconducting 
phase is an island in the x-T plane, completely detached from the T = 0-axis. However, 
this scenario seems hard to realize in real life, as it would require the T = evolution 
curve (the arrow line in Fig. |]) to pass underneath but very close to the PMS phase. 



7 Conclusions 

We have formulated and investigated a simple, spin-Dirac fermion field theory in 2+1 
dimensions capable of explaining, at least qualitatively, a variety of experimental properties 
of the cuprate superconductors and their parent compounds as a function of doping fraction 
and temperature. 

Our model provides a qualitative understanding of insulating, AFM behaviour at low 
doping, of high-T c superconductivity through the Bose-Einstein condensation of spin- 
disorder-bound charge pairs at intermediate doping, and of Fermi-liquid behaviour at 
large doping. In addition we have formulated several predictions, which may be amenable 
to experimental testing. In particular, we recall the possibility of reentrance (Sect. |6.2j ) 
and the statement that superconductivity is unlikely to occur in materials with S > 1/2 
(Sect. |6lD. 

The model is constructed as an effective theory, using as ingredients only points 2,4 
and 7 of the list of important experimental properties given in Sect. |l[ As output, the 
model gives a reasonable explanation for the other points 1,3,5 and 6. We feel this is 
because we have been able to identify the essential degrees of freedom, and to pursue the 
consequences of the relevant symmetries in classifying their possible interactions. 

If we succeed in computing, by means of Monte Carlo simulations, the masses and 
the binding energies present of the various states in the model, we will also be able to 
study crossover from "BE-like" to "BCS-like" behaviour. From the point of view of our 
model, this crossover occurs when going from a situation (i.e., a (k,y) value) in which 
Tbe < Tsg, to a situation in which masses and binding energies are such that Tgc is 
smaller than a would-be Tbe- In the latter case, the quantum liquid condensation would 
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occur simultaneously with the pair formation. 

The answer to these questions, and to many others, will require much more future 
work. However, let us conclude this paper by trying to anticipate our answers to the 
questions listed in Ref. || and quoted in Sect. ||: 

• What is the physical origin of the anomalous normal state? A dynamical, 
antiferromagnetically interacting spin background, strongly coupled to fermions (the 
heavy fermions of the type (r • cf>) if) in our model). 

• What characterizes this anomalous normal state? The presence of heavy 
single fermionic charges ((r • 4>) ip) and light, bosonic charge-2 bound states. 

• What is the mechanism of high-T c superconductivity? Bose-Einstein con- 
densation of these dynamically formed, stable bosonic charge-2 pairs. 

• What is the pairing state? A bound state of heavy fermions, bound by spin-waves 
in a disordered phase: a PMS-pair. 
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